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I .  INTRODUCTION 


« 


t 


1  2 

In  1973,  McClellan  and  Parks  *  published  a  listing  of  a  computer 
program  for  the  design  of  finite-duration  impulse-response  digital 
filters.  This  program  was  unique  in  that  the  authors  had  developed  a 
unified  theory  for  the  design  of  the  four  types  of  filters:  bandpass, 
bandstop,  Hilbert-transform,  and  differentiation.  The  resulting  soft¬ 
ware  is  one  of  the  most  flexible  digital  filter  design  programs  available. 

In  the  analysis  of  ballistic  data  which  has  been  converted  from  an 
analog  voltage  record  to  a  digital  time  series,  it  is  generally  desirable 
to  be  able  to  isolate  various  signal  components  for  individual  study. 

Such  components  are  usually  identifiable  by  frequency  content,  and  as  a 
consequence,  are  ideally  suited  for  separation  or  removal  by  digital 
filtering  techniques.  This  report  describes  the  adaptation  of  the 
filter  design  program  to  the  CDC  computer  at  BRL,  the  modification  of  a 
portion  of  this  program  into  a  subroutine,  the  development  of  subroutines 
to  specify  filter  design  parameters,  and  a  convolution  subroutine  for 
filters  of  odd  length.  A  complete  description  of  design  considerations 
for  digital  filters  and  definitions  of  related  terms  is  beyond  the  scope 
of  this  report.  Any  of  the  cited  references  will  provide  the  necessary 
information.  This  report  does  provide  sufficient  design  information  to 
allow  the  reader  to  implement  digital  filters;  a  subsequent  BRL  Technical 
Report  will  cover  in  greater  detail  specific  application  techniques. 

II.  A  DESCRIPTION  OF  THE  DIGITAL  FILTER  DESIGN  PROGRAM 

Only  minor  changes  were  made  to  the  program  as  it  appeared  in 
reference  2.  The  program  statement  added  was 

PROGRAM  DESIGN (INPUT, OUTPUT, TAPE5=INPUT,TAPE6=0UTPUT,TAPE7) 

In  the  original  program,  when  the  value  of  the  variable  JPUNCH  was  input 
as  1,  values  of  the  filter  coefficients  were  output  to  punched  cards. 

In  the  program,  as  it  exists  on  the  CDC  computer,  TAPE7  may  be  specified 
in  the  jobstream  to  be  any  suitable  device  or  file.  The  constants  PI  and 
PI2  (tt  and  2tt,  respectively)  were  extended  to  the  full  double  precision 
word  length  for  the  CDC.  The  free-field  input  form  of  the  original 
program  was  replaced  with  formatted  input.  Finally,  a  test  for  end-of- 
file  on  input  was  added  to  allow  for  multiple  designs  per  computer  run. 

A  listing  of  program  DESIGN  is  in  Appendix  A. 

^J.H.  McClellan,  T.W.  Parks,  "A  Unified  Approach  to  the  Design  of 
Optimum  FIR  Linear-Phase  Digital  Filters, M  IEEE  Trans.  Circuit  Theory, 
CT-20 (6) ,  697-701  (1973). 

2 

J.H.  McClellan,  T.W.  Parks,  L.R.  Rabiner,  "A  Computer  Program  for 
Designing  Optimum  FIR  Linear  Phase  Digital  Filters,  nIEEE  Trans  Audio 
Electroacoustics,  AU-21(6),  506-526  (1973). 
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III.  A  DESCRIPTION  OF  THE  DIGITAL  FILTER  DESIGN  SUBROUTINE 


For  most  applications  to  analysis  of  ballistic  data,  optimum  digital 
filter  design  specifications  are  the  result  of  a  systematic  trial-and- 
error  investigation.  Frequently,  the  design  specifications  change  from 
one  data  event  to  the  next  because  certain  aspects  of  the  experiment  were 
non-repeatable.  In  view  of  these  factors,  it  seemed  appropriate  to 
formulate  a  filter  design  subroutine  for  use  in  interactive  analysis 
computer  programs,  thereby  permitting  tailoring  of  the  filter  design 
on  a  round-by-round  basis. 

Subroutine  FILTER,  a  listing  of  which  appears  in  Appendix  B,  is 
extracted  from  program  DESIGN.  It  will  design  bandpass  filters  of  up  to 
10  bands,  but  will  not  design  Hilbert  Transformers  or  differentiators. 

The  grid  density  (LGRID)  has  been  fixed  at  16,  but  the  subroutine  other¬ 
wise  retains  the  full  flexibility  of  program  DESIGN.  All  variable  names 
used  in  program  DESIGN  are  also  retained. 

The  subroutine  statement  is 

SUBROUTINE  FILTER(NFILT,NBANDS, EDGE,FX,WTX, IPRINT,H) ,  where  NFILT  is 
the  filter  length;  NBANDS  is  the  number  of  pass/stop  bands;  EDGE  is  an 
array  containing  the  band  edges,  expressed  as  fractions  of  the  sampling 
frequency;  FX  is  an  array  containing  the  desired  filter  shape,  (1.  in 
the  pass  bands  and  0.  in  the  stop  bands);  WTX  is  an  array  containing  the 
desired  relative  weighting  in  each  band;  IPRINT  is  a  control  variable 
for  printing  the  coefficients  (0-print  coefficients,  l-don't  print 
coefficients) ;  and  H  is  the  array  containing  the  filter  coefficients  on 
output.  The  variables  NFILT,  NBANDS,  AND  IPRINT  are  integers;  the  arrays 
EDGE, FX, WTX,  and  H  are  real,  dimensioned  2*NBANDS,  NBANDS,  AND  (NFILT+l)/2, 
respectively.  If  NFILT  is  even,  then  the  H  array  is  dimensioned  NFILT/2. 

IV.  CONSIDERATIONS  IN  THE  USE  OF  THE  DESIGN  SOFTWARE 

For  the  purposes  of  this  discussion,  assume  that  the  data  sequence 
x.  n  consists  of  points  equally  spaced  in  time;  in  particular.  At  will 
1  i=l 

denote  the  time  between  two  consecutive  samples.  The  sampling  frequency, 
fs,  is  therefore  1/At,  and  the  bandwidth  of  the  data  is  .5fs.  The  band¬ 
width  of  the  data  represents  the  highest  unaliased  frequency  present  in 
the  data,  provided  due  care  has  been  given  to  the  sampling  process. 

The  essence  of  the  design  algorithm  is  to  approximate  the  desired 
filter  response  function  on  the  frequency- amplitude  plane  from  -.5fs  to 
+.5fs  on  the  frequency  axis.  The  coefficients  are  designed  in  a  normalized 
form  on  the  interval  [-.5,  .5].  Moreover,  the  frequency  response  has 
either  odd  or  even  symmetry  about  the  origin  on  the  frequency  axis,  so 
that  the  design  problem  is  completely  determined  by  specifying  the 
desired  response  on  the  normalized  frequency  interval  [0.,  .5]. 

In  Figure  1,  below,  is  shown  the  frequency  response  of  a  typical  low 
pass  filter.  This  is  a  two  band  filter:  it  has  a  pass  band  and  a  stop 
band . 
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Figure  1.  Frequency  response  of  a  low  pass  filter 


The  pass  band  is  from  0.  to  .2,  or  to  40%  of  the  bandwidth.  The  frequency 
.2f$  is  termed  the  cutoff  frequency  of  the  filter.  It  is  a  "pass"  band 
since  frequencies  in  this  band  are  "passed"  unaltered  (i.e.  are  multiplied 
by  1).  The  stop  band  is  from  .3  to  .5;  frequencies  in  this  band  are 
"stopped"  (i.e.  multiplied  by  0). 

The  frequency  band  from  .2  to  .3  is  termed  the  transition  band. 
Selection  of  the  width  of  this  transition  band  is  somewhat  critical  in 
the  design  of  a  digital  filter,  for  the  following  reason:  as  the  transi¬ 
tion  band  narrows,  the  slope  of  the  frequency  response  (_i.e.  the  filter 
roll-off)  increases.  As  this  slope  increases,  the  design  algorithm 
compensates  by  increasing  the  deviation  from  the  desired  response  in  the 
pass  and  stop  bands.  This  deviation  is  called  the  "ripple",  and  results 
in  increases  and  decreases  of  amplitude  in  the  signal  at  those  particular 
frequencies.  An  example  of  a  filter  designed  with  too  narrow  a  transi¬ 
tion  band  is  shown  in  Figure  2. 

In  any  application  software,  it  is  advisable  to  have  the  capability 
of  viewing  the  frequency  response  of  the  filter  prior  to  its  application, 
in  order  to  be  certain  of  its  characteristics.  The  design  program,  as 
a  part  of  its  printed  output,  lists  the  normalized  frequencies  at  which 
the  maximum  and  minimum  amplitudes  of  the  ripple  occur.  Also  listed  are 
the  deviations  from  the  desired  design,  which  provide  a  measure  of  the 
amplitude  error  to  be  expected  as  a  result  of  applying  the  filter  to  the 
data.  (See  Appendix  C). 

Referring  to  the  example  of  Figure  1,  the  input  variables  to  design 
this  filter  were  assigned  the  following  values: 
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or 


Figure  2.  Frequency  response  of  a  low  pass  filter  with  narrow 
transition  band 


NFILT  =  33 

NBANDS  =  2 
EDGE ( 1 )  =  0 
EDGE(2)  =  .2 
EDGE(3)  =  .3 
EDGE ( 1 )  =  .5 
FX(1)  =  1. 

FX ( 2 )  =  0. 

WTX(l)  =  10. 

WTX(2)  =  100. 

While  NFILT  is  specified  as  33,  only  17  distinct  coefficients  are  returned, 
since  the  design  is  symmetric  about  0.  The  sample  output  in  Appendix  C 
indicates  the  ordering  of  the  17  coefficients,  although  this  is  not  the 
fi  1  ter  of  Fi  gure  1 . 

As  can  be  seen  in  this  example,  the  EDGE  array  specifies  the 
normalized  band  edges.  The  FX  array  specifies  the  desired  amplitude 
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A 

v/ 


response,  which  is  usually  (but  not  necessarily)  1  in  the  pass  bands  and 
0  in  the  stop  bands.  The  WTX  array  specifies  a  relative  scaling  of  the 
magnitude  of  the  deviation  between  the  pass  band  and  the  stop  band.  In 
this  example,  the  design  algorithm  allows  10  times  less  deviation  in  the 
stop  band  than  in  the  pass  band.  This  relative  weighting  may  be  adjusted 
arbitrarily  to  suit  a  particular  need.  For  example,  by  relaxing  the 
pass  band  weighting,  say  WTX(1)-1.  ,  one  could  design  a  filter  with  a  more 
narrow  transition  band. 


For  additional  information  concerning  the  design  of  digital  filters, 
the  reader  is  referred  to  references  1,2,  and  3. 

V.  IMPLEMENTATION  OF  DIGITAL  FILTERS 


A  digital  filter  is  applied  to  a  data  sequence  by  convoluting  the 
filter  weights,  or  coefficients,  with  the  data  points.  Specifically,  if 


Mi.,  • 


is  a  sequence  of  data  points  equally  spaced  in  time,  and  if 


ML 


are  the  filter  coefficients,  where  N_<n,  then  the  filtered  data  sequence 

iy.V1  has  values  given  by 
\  1/i=N 


yi  = 


X 

k=l 


\Xi+l-k 


(1) 


One  notes  that  if  i<N-l,  then  i+l-k<N-k,  so  that  some  subscripts  of  x  in 
the  summation  may  have  values  less  than  1;  we  have  no  corresponding  x 
values.  There  are  two  choices:  either  start  the  convolution  process  at 
i-N,  or  start  at  i=l  and  modify  k  to  avoid  subscripts  of  x  less  than  1 
until  we  get  to  the  Nth  point.  In  the  first  case,  N-l  data  points  at  the 
beginning  are  unused,  and  the  output  sequence  starts  at  i-N.  In  the  second 
case,  the  first  N-l  output  points  have  not  been  transformed  by  the  same  set 
of  filter  coefficients  as  have  the  rest  of  the  data;  the  first  N-l  points 
are  of  questionable  value.  It  will  be  noted  that  the  same  problem  occurs 
for  i>n-N.  In  what  follows,  we  will  discuss  a  method  to  avoid  these  diffi¬ 
culties.  In  particular,  it  will  be  shown  that  an  n-point  input  sequence 
can  be  modified  to  provide  n  useful  output  points. 


For  nonreal-time  applications,  i.e.:  for  post-processing  of  data, 
one  is  in  the  admirable  position  of  knowing  in  advance  what  is  going  to 
occur.  That  is,  the  convolution  process  can  be  numerically  manipulated 
so  as  to  provide  one  output  point  corresponding  to  each  input  point, 
with  no  lag.  (Only  filters  of  finite  odd  length,  say  N  =  2M  +  1,  M  = 
1,2,..., 63,  will  be  discussed  here.)  This  is  accomplished  simply  by 
moving  the  filter  coefficients  M  indeces  in  Eq.  (1),  so  that 


yi = 


2M+1 

X 

k=l 


hkXi-M-l+k 


(2) 


Eq.  (2)  implies  that  the  filter  coefficients  are  centered  at  the  ith  data 
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f  I  1 

point.  If  the  coefficients  are  re-indexed  as  ,  then  Eq .  (2)  is 

more  simply  written  as 


k=-M 


where  h, 

k 


^k+Mfl 


(3) 


Now,  for  the  values  i=l, 2 , . . . ,M,n~M+l ,n-M+2, . . . ,n-l,n,  Eq .  (3)  still  has 
some  values  of  i-k  for  which  there  is  no  corresponding  x. 

It  is  necessary  to  provide  M  values  at  the  beginning  and  M  values  at  the 

end  of  the  sequence  /x_^n  •  This  can  be  done  with  a  minimum  of  frequency 

distortion  by  using  an  oid  reflection  of  the  first  M  and  last  M  points.  Speci¬ 
fically,  for  i-k<l,  define  x^k  by 


x 


i-k 


Similary , 


=  2xrxk-±+2  • 

for  i-k>n,  define  x^_k  by 


(4) 


x 


i-k 


2x  -x 
n 


k-i 


(5) 


Graphically,  Eq.  (4)  reflects  X2,X3’**#XM  about  the  vertical  line 
through  x1  and  then  about  the  horizontal  line  through  xx.  The  points 

x  x  , . . .x  are  reflected  in  a  like  manner  about  x  ,  as  shown 

n-M+1 5  n-M+2 5  n-1  n 

in  Figure  3. 


This  reflection  process  can  be  trivially  incorporated  into  the 
convolution  algorithm,  as  will  be  explained  below. 

The  types  of  digital  filters  being  considered  here  have  an  additional 
property  which  simplifies  the  convolution  process:  they  are  of  either 
even  or  odd  symmetry  about  their  midpoint.  That  is, 

h  =  ±  h  ,  k=l, 2, .  .  .M  (6) 

k  -k 
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Figure  3.  Graphical  construction  of  new  end  points  for  x_^ 


As  a  consequence,  Eq .  (3)  may  be  written  as 

M 

y .  =  h  x .  +  2  Ir  (x .  ±  x  .  , .  ) 

o  1  k  l-k  i+k 

k=l 


(7) 


Whereas  Eq.  (3)  requires  2M  +  1  multiplications  and  2M  additions  to 
implement,  Eq .  (7)  requires  only  M  4-  1  multiplications  and  M  additions. 

Utilizing  Eqs.  (4),  (5),  and  (7),  the  following  are  three  examples 
of  convolution  subroutines.  The  first,  in  FORTRAN,  is  in  use  on  the  BRL 
CDC  system.  The  second  is  in  standard  BASIC.  The  third,  in  an  enhanced 
BASIC,  is  in  use  on  several  BRL  systems.  In  each  case,  x  is  the  input/ 
output  array  of  length  N.  The  K=M+1  filter  coefficients  are  stored  in 
the  array  H. 

Example  1:  FORTRAN  Convolution  Subroutine 

SUBROUTINE  CONVOL  (H,K,X,N) 

DIMENSION  H(K) ,X(N) , S (127) , SAVE (63) 

M  =  K-l 
L  =  K  +  M 
IF (L. GT .127)  STOP 
DO  5  1=1, M 

S(I)  =  2.*X(1)  -  X(K  -I) 
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S(L+1-I)=X(I) 

SAVE (I)  =  2.*X(N)  -X(N-I) 

5  CONTINUE 
S(K) =X (1) 

LAST=N-M 
DO  20  1=1, N 
X(I)  =  0. 

DO  10  J=1,M 

X(I)=X(I)+H(J)*(S(J)+S(L+1-J) 

10  CONTINUE 

X(I)=X(I)+H(K)*S(K) 

DO  15  J=2,  L 
S(J-1)=S(J) 

15  CONTINUE 

lF(I.LE.LAST)  S(L)=X(I+K) 
lF(I.GT.LAST)  S(L)=SAVE(I-LAST) 
20  CONTINUE 
RETURN 
END 


Example  2:  BASIC  Convolution  Subroutine 

100  SUBROUTINE  Convolution  (H,K,X,N) 

110  DIM  H(K) ,X(N) ,S(127) ,SAVE(63) 

120  M=K-1 

130  L=K+M 

140  IF  L>127  THEN  STOP 

150  FOR  1=1  TO  M 

160  S(I)  =  2.*X(1)-X(K-I) 

170  S(L+1-I)=X(I) 

180  SAVE(I)=2.*X (N)-X(N-I) 

190  NEXT  I 
200  S (K) =X(1) 

210  Last-N-M 
220  FOR  1=1  TO  N 
230  X(I)=0 

240  FOR  J=1  TO  M 

250  X(I)=X(I)+H(J)* (S(J)+S (L+l-J) ) 

260  NEXT  J 

270  X(I)=X(I)+H(K)*S(K) 

280  FOR  J=2  TO  L 
290  S(J-1=S(J) 

300  NEXT  J 

310  IF  I<=Last  THEN  S(L)-X(I+K) 

320  If  I>  Last  THEN  S (L)=SAVE(I-Last) 
330  NEXT  I 
340  RETURN 
350  SUBEND 
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In  the  third  example,  use  is  made  of  several  matrix  operations  available 
in  enhanced  BASIC.  The  function  DOT  returns  the  dot  product  of  the  two 
input  arrays.  The  function  MAT  REORDER  rearranges  the  elements  of  one 
array  according  to  the  index  order  specified  by  another.  In  this  example, 
the  array  B  has  the  values  2,3,4, ... ,L,1,  where  L=2M+1,  the  filter  length. 
Implementation  of  the  routine  in  example  3  represents  a  decrease  in 
execution  time  by  a  factor  of  15  over  the  routine  in  example  2. 

Example  3:  BASIC  Matrix  convolution  Subroutine 

100  SUBROUTINE  Convolution (H,L,X,N, B) 

110  DIM  H (L) ,X(N) , S (127)  ,  Save(63) ,B (L) 

120  IF  L>127  THEN  STOP 
130  REDIM  S(L) 

140  M=  INT(L/2) 

150  FOR  1=1  TO  M 

160  S(I)=2.*X(1)-X(M+1-I) 

170  S(L+1-I)=X(I) 

180  Save(I)=2.*X(N)-X(N-I) 

190  NEXT  I 
200  S(M+1)=X(1) 

210  Last=N-M 
220  FOR  1=1  TO  N 
230  X(I)=DOT (H, S) 

240  MAT  REORDER  S  BY  B 

250  IF  K=Last  THEN  S  (L)=X(I+M+1) 

260  IF  I>Last  THEN  S (L)=SAVE (I-Last) 

270  NEXT  I 
280  RETURN 
290  SUBEND 


VI.  CONCLUSIONS 

Digital  filters  have  a  wide  range  of  application  for  numerical 
analysis  of  time-series  data.  The  filter  design  program  presented  here 
has  been  found  to  be  one  of  the  most  versatile  available.  The  reflection 
principle  described  in  this  report  seems  to  introduce  the  least  additional 
frequency  content  into  the  data  of  any  of  the  methods  available.  This 
same  technique  has  been  used  to  produce  periodic  continuation  of  essentially 
transient  phenomena,  facilitating  the  use  of  numerical  filters  in  their 
analysis . 

In  a  forthcoming  BRL  Technical  Report,  the  author  will  discuss 
specific  techniques  for  the  application  of  digital  filters  to  the  analysis 
of  ballistic  data.  The  report  will  also  develop  in  greater  detail  the 
theory  and  applicability  of  digital  filters  to  analysis  of  time  series. 

VII .  SUMMARY 

An  open  literature  FORTRAN  computer  program  for  the  design  of  finite 
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impulse-response  digital  filters  has  been  implemented  on  the  BRL  CYBER 
system.  Algorithms  have  been  developed  and  coded  for  the  convolution  of 
digital  filters  with  time  series  data.  These  algorithms  include  a  method 
for  the  removal  of  the  filter  delay,  as  well  as  elimination  of  the  loss  of 
data  at  the  beginning  and  end  of  the  particular  data  set  being  filtered. 
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APPENDIX  C 

SAMPLE  OUTPUT  FROM  PROGRAM  DESIGN 
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USER  EVALUATION  OF  REPORT 


Please  take  a  few  minutes  to  answer  the  questions  below;  tear  out 
this  sheet,  fold  as  indicated,  staple  or  tape  closed,  and  place 
in  the  mail.  Your  comments  will  provide  us  with  information  for 
improving  future  reports. 

1 .  BRL  Report  Number _ _ 

2.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related 
project,  or  other  area  of  interest  for  which  report  will  be  used.) 


3.  How,  specifically,  is  the  report  being  used?  (Information 
source,  design  data  or  procedure,  management  procedure,  source  of 
ideas,  etc.) _ 


4.  Has  the  information  in  this  report  led  to  any  quantitative 
savings  as  far  as  man-hours/contract  dollars  saved,  operating  costs 
avoided,  efficiencies  achieved,  etc.?  If  so,  please  elaborate. 


5.  General  Comments  (Indicate  what  you  think  should  be  changed  to 
make  this  report  and  future  reports  of  this  type  more  responsive 
to  your  needs,  more  usable,  improve  readability,  etc.) _ 


6.  If  you  would  like  to  be  contacted  by  the  personnel  who  prepared 
this  report  to  raise  specific  questions  or  discuss  the  topic, 
please  fill  in  the  following  information. 

Name : 


Telephone  Number: 
Organization  Address: 


